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Abstract. We derive steady equilibria for lateral downslope moisture flow in an ide- 
alized thin closed layer as a solution to the ID Richards' Equation. The equilibria are 
determined by two free parameters: the downslope flux and a boundary condition. So- 
lutions exhibit a constant downslope flow speed and moisture content for the constant 
equilibrium flux, which is the product of the two. However where an isolated zone of fixed 
saturation degree exists representing a boundary condition, the flow speed immediately 
upslope is reduced and the moisture content correspondingly increased to preserve the 
constant equilibrium flux. The capillary head jump at the saturated zone produces a block- 
age that gives a high moisture content back upslope through a pooling distance deter- 
mined by the equilibrium condition that the downslope flux is constant. In our numer- 
ical integrations, the vertically projected pooling height is more than 10 km for a fully 
saturated zone in mixed silty or clay soils, but decreases by about an order of magni- 
tude with every 10% decrease in the boundary-zone saturation degree. The drying of down- 
hill saturated zones with the increased speed of mountain moisture outflow and corre- 
sponding decreased mountain moisture content gives a viable explanation for the mys- 
terious 69% unaccounted drop seen in the spring outflow in the La Luz / Fresnal Wa- 
tershed at Alamogordo's upstream spring-box diversions in the semiarid southeastern New 
Mexico USA. 

[agu index terms: hydrology (1800), groundwater hydrology (1829), groundwater trans- 
port (1832), streamflow (1860), soil moisture (1866), vadose zone (1875)] 



1. Introduction 

Largely with the hope of increasing the captured outflow 
of the La Luz / Fresnal Canyon stream system, the Town of 
Alamogordo New Mexico substantially repaired and moved 
its collection system between 1982-1997, to divert nearly 
all of the streamflow at or just below the mountain spring 
sources. The locale is shown on the topo map in Figure 1 
with the main springs indicated. By transporting the water 
in pipes through the lower canyons of the western Sacra- 
mento Mountains they minimize evapotranspiration losses 
and avoid possible losses due to unappropriated human us- 
age, which is believed to have been a problem historically. 

The 100-year-span USGS-monitored baseflow at the foot 
of the mountains near the village of La Luz was 7.8 million 
gallons per day (Mgal/d) (342 1/s), shown as a dashed line in 
Figure 2a (USGS reports summarized in Table E-6 of South 
Central Mountain RC&D Council et al. 2002), but in the 
mid 1980s the streamflow ranged up to about 11.5 Mgal/d 
(504 1/s) in the median of the daily averages selected from 
two-month periods when the flows were consistently high 
(plus signs, from New Mexico Water Resources Data, USCS 
yearly reports). The La Luz / Fresnal streamflow was quite 
variable since the USGS monitoring resumed after 1980 due 
to the intermittent operation of Alamogordo's diversions, 
but during 1985-87 the La Luz spring capture boxes were 
not operated and the streams allowed to flow fully during 
the transition to the new system. A consistently high flow 
during the transition period indicates that Alamogordo was 
not diverting the water elsewhere in the system and allows 
us to obtain a best minimum estimate for the full natural 
streamflow. The empty boxes in Figure 2a show the av- 
erage metered pipe flow made during periods when some 
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streamflow was allowed and give a lower bound on the ac- 
tual system outflow; the filled boxes are for periods when 
the entire flow was piped and streams dryed, and represent 
essentially the full system outflow (New Mexico Office of 
the State Engineer, unpublished data, 2006; Town of Alam- 
ogordo, unpublished data, 2007) . The data point rms accu- 
racy is limited by the inherent device precision, which for 
Alamogordo's pipe-flow measurements is 2%, and for the 
USGS streamflow measurements is estimated to be 5%, the 
minimum noise for the daily averages obtained with their 
system. 

However after the upstream diversion relocation was com- 
plete, rather than increase, the total outflow of the system 
dropped to less than 2.5 Mgal/d (110 1/s) by 2002. The drop 
represents a huge decrease from about 11.5 Mgal/d (504 1/s) 
to less than 2.5 Mgal/d (110 1/s) or a loss of about 78%! Dur- 
ing the same 1986-2002 period, the average precipitation in 
the Cloudcroft and mountain source area decreased by about 
25%, as shown by the smoothed curve in Figure 2c ( West- 
ern Regional Climate Center 2006). Similar Tularosa Creek 
showed a larger 30% drop during the period suggesting a 
streamflow overresponse to the precipitation change and im- 
plying a reduced projected La Luz / Fresnal flow for the dry 
period of about 8 Mgal/d (351 1/s). Thus without other sig- 
nificant appropriations or diversions, we estimate a drop in 
the total system outflow from about 8 Mgal/d (351 1/s) to 
less than 2.5 Mgal/d (110 1/s) meaning that about 69% of 
the normal system outflow cannot be accounted for since the 
mid-1980's. The unaccounted drop is probably conservative 
since downstream losses are not considered in comparing up- 
stream spring outflow to downstream flow. However small 
secondary sources are also known to exist, which act oppo- 
sitely, adding to the downstream flow. 

As a baseline, the USGS daily record from Bent NM 
shows that Tularosa Creek about 10 miles (16 km) to the 
north, which has remained intact, has changed consistently 
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Creek system, and apparently nearly uniform loss in all of 
the springs after 1997 when Alamogordo started taking all of 
the streamflow, suggests a hydrological phenomenon. Rick 
Warnock, president of the Sacramento Mountain Watershed 
Restoration Corporation, has maintained that the lower- 
canyon drying must be the main cause of the mountain wa- 
ter depletion (Warnock July 11, 2006). In one study in an- 
other system, diurnal and bi-annual moisture variations are 
found to be strongly linked to changing downstream river 
flow (Waichler and Yabusaki 2005). 

Such a large anomalous event gives unique opportunity 
for a critical hydrological test: Is it possible that down- 
stream drying could cause upslope moisture depletion? A 
downstream influence on upslope conditions may be possible 
if unsaturated lateral downslope moisture dynamics plays a 
dominant role in controlling the water distribution through- 
out the system. The relative importance and interrelation of 
saturated/unsaturated dynamics in downslope flow has been 
a matter of study and discussion (Freeze and Witherspoon 
1966; Jackson 1992). Numerous measurements indicate di- 
rect 3D unsaturated near-surface moisture flow, suggesting 
that a significant downslope flow component may remain as 
unsaturated (e.g. Weyman 1973; Nieber and Walter 1981; 
Miyazaki 1988; Stnai and Dirksen 2006). A largely unsatu- 
rated downslope flow may be suggested by the fact that the 



Figure 1. Topo Map of western Sacramento Mountains 
NM USA with contour interval 250 ft (76.2 m) and heavy- 
contour interval 1000 ft (305 m) from 5000 ft (1524 m) 
to 9000 ft (2745 m) above sea level; the total vertical 
range is 4200 ft (1280 m); blue dots are the main springs; 
normal streams are shown as solid blue and normally dry 
streambeds as dashed blue; the black dashed line denotes 
the foot of the mountains; major roads are shown as red. 



with the mountain precipitation with a 3.76 year delayed 
response, evident in comparing the smoothed data lines be- 
tween Figures 2b and 2c (New Mexico Water Resources 
Data, USGS yearly reports). Since 1998, Cloudcroft and its 
surroundings above the La Luz / Fresnal Watershed have 
experienced considerable drying of wells and springs, with a 
forest mortality consistent with a reduced mountain mois- 
ture content, whereas no similar depletions or its effects 
seem evident in Mescalero and the mountains above Tu- 
larosa Creek. Laborcita Creek, which has a much smaller 
normal outflow, was not monitored consistently during the 
period. 

The drop in outflow seen in the La Luz / Fresnal Stream 
System is regionally localized, appears to be historically un- 
precedented, and has been persisting after 2002, so most 
common effects can be ruled out. For example, precipita- 
tion patterns, snowpack variation or storm distribution and 
frequency, have larger spatial scale and are historically re- 
curring. A similarly large drop in outflow was not seen in 
the other major stream systems in the locale, and an un- 
usually large drop did not occur in the La Luz / Fresnal 
Stream System following the 1970s dry period. The loss 
is greater than all the other appropriations in the system 
combined and much larger than the 0.56 Mgal/d (24.5 1/s) 
increase in appropriations for the 1980-2000 period, which 
was for domestic wells mainly located below the Alamogordo 
diversions (from New Mexico Office of the State Engineer, 
WATERS database www.ose.state.nm.us, 2007). Other wa- 
ter users in the system have experienced similar drops in 
their spring flow and drying of their wells during the period 
too. Thus a significant impact due to other human usage 
seems precluded. 

The steady and precipitous drop with approximately the 
same e-folding time as the delayed response in the Tularosa 
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Figure 2. (a) La Luz / Fresnal System partial- (empty 
boxes), and full- diversions (filled boxes), historic aver- 
age baseflow (dashed line), and medians of daily flows 
(plus signs); (b) Tularosa Creek yearly medians (boxes), 
smoothed with a 4-year FWHM Gaussian (line), and his- 
toric average baseflow (dashed); (c) Cloudcroft precipita- 
tion yearly averages and 4-year smoothed (line). 
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Figure 3. Idealized ID friction flow down a uniformly 
filled closed flow layer of constant thickness d and length 
L, inclined at i from the horizontal x, repeated infinitely 
into the plane of the drawing in y; z is the vertical and I 
the downslope coordinate, which ranges from to L. 
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streamflow in these systems is never larger than about 15% 
of the total precipitation inflow. 

We seek guidance on the gross characteristics of large- 
scale mountain moisture storage dynamics outlined here. 
The geological and topographical constraints suggest that 
moisture leaves the mountains traveling through the lower 
canyons mainly as a lateral downslope flow in a thin 
layer being prevented from seeping down very far be- 
cause of underlying impervious layers. The Richards' 
Equation describes general moisture flow in mixed satu- 
rated/unsaturated ground- water systems giving a correct 
upper vadose zone, capillary fringe, water table, and lower 
saturated zone in simulations, though its multidimensional 
form is difficult to model due to its nonlinear properties 
(Fipps and Skaggs 1989; Fiori and Russo 2007). 

One published solution to the Richards' Equation for 
general downslope flow in a thin layer assumes a constant 
downslope flow speed and moisture content (Philip 1991), 
while others show pooling at the bottom of the hillslope 
(Sloan and Moore 1984; Hurley and Pantelis 1985; Stagnitti 
et al. 1986; Steenhuis et al. 1999). As we show in this study, 
flow solutions of both these types are obtained depend- 
ing upon the choice of downhill boundary condition. The 
boundary condition may be affected by streamflow main- 
taining a downhill saturated zone. Where a high down- 
hill saturation degree is maintained, moisture backs upslope 
through a pooling distance that depends upon the soil type. 
Alamogordo's diversions near the canyon spring sources dry 
the lower-canyon stream system, which must dry the lower 
saturated zones and change the downhill boundary condi- 
tion. Such drying leads to a higher-speed lateral downs- 
lope moisture flow into the open alluvial basin below with 
a decreased moisture content back upslope as we describe. 
The reduced moisture content can affect the nearby spring 
outflow as well as possibly the relative mountain moisture 
content above. 

2. ID Steady Moisture Flow 

Friction flow down a uniformly filled closed layer of con- 
stant inclination i and thickness d is an idealization of the 
flow of mountain moisture to a lower basin outlet over a dis- 
tance L, as illustrated in Figure 3. By Darcy's rule, the flux 
(or specific flux) q with units of speed (L/T) is always away 
from the pressure head (Chapter 5, Bear 1988) 



q= -K(S) (V^(S) + Vz 



(1) 



written for isotropic conditions, where the scalar hydraulic 
conductivity K(S) with units of speed (L/T) and the cap- 
illary pressure head tp(S) with units of length (L) are func- 
tions of the saturation (or saturation degree) S and the soil 
properties; Vz represents the gradient of the gravitational 
pressure head. The saturation degree S ranges from to 
1 in proportion to the fractional water filling of the avail- 
able spaces between soil grains, with S — 1 indicating fully 
saturated. The flux q is proportional to the flow speed V 



q = 6V = nSV, 



(2) 



where 8 = nS defines the moisture content (Section 9.4.4 
Bear 1988), and n is the porosity, the fractional volume of 
open space between soil grains. Though in our mathematical 
formulation a spatially varying porosity is allowed, in gen- 
eral discussions concerning the downslope flow speed and 
moisture content, a uniform porosity is assumed for simplic- 
ity. 

The conservation of moisture is represented by the incom- 
pressible form of the continuity equation (Section 9.4, Bear 
1988) 



dS 



(3) 



assuming that the medium is nondeformable or the porosity 
n not temporally varying. We have neglected other terms 
commonly introduced into the continuity equation repre- 
senting precipitation influx, evapotranspiration loss, or seep- 
age out the bottom below the hilltop entry point at £ — 0. 
While such effects might be included, they should not al- 
ter the essential physical properties of the solutions, as we 
discuss further in Section 7. 

Combining (1) and (3) yields the Richards' general flow 
equation (Richards 1931) 



ds ~ 



(k(S)(^VS + Vz)), (4) 



having expanded the gradient of the capillary head with its 
implicit spatial dependence Vip(S(x)) = (dtp / dS)V S (x) . 

We are primarily interested in the lateral downslope com- 
ponent to the flux qe, which satisfies the £ vector-element 
equation of (1) 



= K(S) (s 



d4>_dS_ 
dS d£ 



)• 



(5) 



using dz/dl = — sint, for the hillslope inclination i. 

We adopt the ID approximation as justified and used in 
other similar studies (Sloan and Moore 1984; Hurley and 
Pantelis 1985; Stagnitti et al. 1986; Steenhuis et al. 1999), 
supposing that moisture changes normal to the surface are 
small compared to lateral changes down the flow layer, which 
should be the case when the idealized flow layer of F_igure 3 
is thin. The ID flow vector can be written q = qe\7£, and 
with variations only in the downslope direction £, the flux 
and saturation written qi = qe(£,t) and S = S(£,i). The 
continuity equation (3) in the ID approximation is written 



8S_ _ dq e 

U dt d£ 



(6) 



Steady equilibrium solutions are obtained with dS/dt = 
in (6), which gives a constant flux qi downslope in £. A 
constant saturation down the flow layer S(£) = S for the 
reference saturation degree S gives the constant flux qt and 
represents one steady solution to (5), written 



qt = K(S) sin l. 



(7) 



However general steady equilibria with a varying S(£) and 
nonzero derivative dS/d£ in (5) also exist. Solving for the 
derivative gives the governing differential equation for satu- 
ration degree as a function of position down the flow layer 
S(£) 



dS 



(dAb 

\dS 



1 

\ K(S))' 



(8) 



now using the reference saturation degree S as a problem 
constant to represent the constant flux q e defined by (7). 
Since dz — —d£smi, dependencies on the inclination angle 
simply project, and the equation exhibits universal solutions 
S(z) as functions of the vertical coordinate z measured pos- 
itively upward and independent of the inclination angle i 



dS = (d4 L y 1 (K(S) 
dz \dS) \K(S) 



(9) 



Solutions S(z) to the first-order differential equation are de- 
scribable by one additional free parameter besides the con- 
stant flux qe or reference saturation S, a boundary value 
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Figure 4. Van Genuchten material functions with sat- 
uration degree S: (a) relative hydraulic conductivity 
K(S) / K ea x\ (b) capillary head tp(S) (m); and (c) deriva- 
tive dtp/dS (m), for: Pachappa fine sandy clay (77 = 1.860, 
tpo = 0.80 m, dotted), Pachappa loam (77 = 2.195, tpo = 
1.65 m, long dash), hypothetical mix 1 (77 = 2.7, "00 = 
1.5 m, solid), mix 2 (77 = 3., Vo = 1-3 m, dashed), mix 3 
(77 = 3.5, t/>o = 1.1 m dash dot), and Amarillo silty clay 
loam (77 = 4.510, tpo = 1.015 m, long-dash dot dot) all 
with S T = 0.1. 



Solutions to (9) are independent of the saturated hy- 
draulic conductivity K sat , which divides out in the ratio 
K(S)/K(S). Such independence gives the equilibria cer- 
tain robustness as K Bat is the least certain material constant 
due to the inhomogeneous character of real ground material, 
which can allow some relatively high-speed horizontal flow 
in cracks and open zones (Sigda and Wilson 2003). It ranged 
from about 0.01 to 3 m/d in field tests (from the averages 
in Table 1 of Xiang et al. 1997), but to much larger values 
> 50 m/d with a median of 7 m/d for a sample of wells from 
draw-down measurements made in the Tularosa Basin (Orr 
and Myers 1986). We suppose K e3 x = 7 m/d for estimating 
downslope flow speeds. The choice gives a saturated flow 
speed Vsat = ^sat sint/n = 5.0 m/d, combining (7) and (2) 
for S — 1, with l — 7° as the half-mountain-height inclina- 
tion, and with the choice n = 0.17 for the porosity of poorly 
sorted soils (Section 2.5.2, Bear 1988). For an approximate 7 
km mid-mountain characteristic downslope spatial scale for 
upper La Luz Creek from Figure 1, we obtain a saturated- 
flow timescale of 3.8 years for the system, consistent with 
the e-folding time for the drop after 1997 seen in Figure 2 
or the delayed response for Tularosa Creek with changing 
precipitation. 




S(£) at some downslope position I, which may represent the 
saturation in a fixed constant boundary zone. The boundary 
condition only makes physical sense as a downhill condition 
given the form of the solutions found in Section 3, as we 
discuss further in Section 6. 

We adopt the van Genuchten material functions for K(S) 
and ip(S) (van Genuchten 1980; van Genuchten and Neilsen 
1985), written in terms of the effective saturation degree 
S e (S) = (S — Si)/(1 — Si), for Sr a small retained moisture 
saturation degree 



k(s) = /Gats c (s)s (!-(!- Sc(s^yy 



i- v 



tp(S) = tpo 



(10) 



(11) 



,Sc(S)- J 

where K s3 x, the saturated hydraulic conductivity, 77, and tpo 
are material constants. The retained moisture saturation 
degree S T is a small zero offset for the saturation curves 
that depends upon the soil type. For our calculations, wc 
suppose uniform soil properties with constant material pa- 
rameters -Ksat, 77, tpo, and S T , and use S T = 0.1 . 

The Van Genuchten formulation is empirical, but is a 
good approximation for soil measurements as described by 
Assouline and collaborators in their fully developed theoret- 
ical formulation of the problem. Assouline and Tartakovsky 
[2001] give best fits for the van Genuchten coefficients for 
different materials in their Table 2. Figure 4 illustrates the 
van Genuchten material functions for various soil samples. 
The singular behavior in the derivative dtp/dS approach- 
ing full saturation S — > 1 leads to a very small saturation 
change with height near full saturation in (9) and solutions 
with large pooling heights. The La Luz / Fresnal Watershed 
consists of mixed soils and should correspond to an interme- 
diate silty soil between coarse sand 77 < 2.2 and dense clay 
77 > 3.5, taken between Pachappa Loam and hypothetical 
mix 2. 
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Figure 5. Saturation degree S as a function of vertical 
drop £sint = —z (km) down the flow layer from numeri- 
cal integrations of the equilibrium equation (9) for a flux 
corresponding to the reference saturation degree S = 0.30 
with different downhill boundary-zone saturation degrees 
(labeled at right). 
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3. Numerical Solutions 

Figure 5 presents equilibrium solutions to (9) formed by 
fourth-order Runge-Kutta numerical integration in the +z 
upslope direction from a given fixed downhill saturation de- 
gree S at z — 20 km, as labeled on the right of each example, 
for the reference saturation degree S = 0.30 and materials 
defined by the van Genuchten coefficient combinations 77 and 
^0 given with Figure 4. The integration step of 1 m was cho- 
sen as the largest step size that always produces stable and 
consistent results. The models illustrate the tendency for 
moisture to backup upslope at near the saturation degree 
of the fixed downhill zone through a characteristic pooling 
height determined by the boundary value and material type. 

Figure 6 shows the pooling height as a function of the 
boundary saturation degree, which is defined by locating the 
intercept of the saturation curve above the lower boundary 
with the upper-hillslope reference saturation line in numer- 
ical solutions as from Figure 5. The pooling height shows 
about a factor of 10 drop with every 10% decrease in the 
downhill saturation degree for any soil type near high sat- 
uration, though clays and silts require less downhill satura- 
tion for a given pooling height. The pooling height is more 
than 10 km behind a downhill saturated zone in all but the 
sandiest soil type, the Pachappa fine sandy clay. For aver- 
age mix 1 or sandier, it decreases to less than 2 km with a 
downhill saturation degree of S = 0.95, and less than a km 
with S < 0.9. Though not shown, pooling heights^ decrease 
a little with smaller reference saturation degree S and are 
not very sensitive to the retained moisture saturation degree 

Sr. 



downslope flux. The flux qe is defined by the product of 
the saturation degree and flow speed in (2), or for the ID 
solutions 



nS(l)Ve{£) = nSVt. 



(12) 



The ID equilibrium solutions from (9) determine the downs- 
lope saturation function S(£) for a given flux qe or refer- 
ence saturation degree S. Given a saturation function S(£), 
we obtain immediately the general product containing the 
downslope flow speed nVe(£) using (12), for what may be 
possibly a downslope varying porosity n(£). Since K(S) is a 
monotonically increasing function of S for all material types 
from Figure 4a, a constant reference saturation degree S is 
uniquely defined by the constant flux qe in (7), so both the 
constant product with the flow speed nVe and constant sat- 
uration degree S for the upper hillslope are determined by 
the given flux qe using (7) and (12). 

The reference saturation degree S must always increase 
with increasing flux qt from (7). However the upper-hillslope 
flow speed in nVe may either increase or decrease with in- 
creasing qe from (12). The rate of change dnVe/dqe can be 
determined by taking the derivative d/dqe in (12), which 
gives 



dqe 



nVt — . 
dqe 



(13) 



then substituting the derivative dS /dqe from (7) and re- 
using (12) and (7) gives 



4. Upper-Hillslope Flow Conditions 

The Richards' Equation solution properties, the varia- 
tions of saturation degree and flow speed, are defined by the 
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Figure 6. Pooling height (km) plotted on a log scale as a 
function of the saturation degree S in the downhill zone, 
synthesizing the numerical experiments for the reference 
saturation degree S = 0.30. 




Figure 7. Sensitivity of upper-hillslope flow speed Ve to 
changing flux qe, dnVe/dqe as a function of the reference 
saturation S for the van Genuchten materials defined in 
Figure 4. To the right of the intercept with the zero 
derivative line (horizontal dashed line), the flow speed Ve 
decreases with increasing flux qe- 



K(S) ( dK(S) \ lN 
S \ dS ) , 



(14) 



The discriminator dnVe/dqe is determined by the hydraulic 
conductivity function K(S) from (10) alone and is plotted in 
Figure 7 for all of the van Genuchten material types defined 
in Figure 4. Where dnVe/dqe is negative, the flow speed in 
the upper hillslope decreases with increasing flux. 

5. Flux Estimate 

Our choice of reference saturation degree S = 0.30 for 
the La Luz / Fresnal Watershed is reasonable considering 
the drop in outflow seen at the upstream spring boxes, but 
seems justified too given estimates of the precipitation in- 
flow. The constant flux down the layer is defined by (7). 
For S = 0.30, l = 7°, and the choice K sa , t = 7 m/d, with 
Pachappa loam as defined in Figure 4, we obtain the flux 
qe = 0.249; with mix 1, qe = 0.326 m/d; and with mix 2, 
qe = 0.354 m/d. 

The specific flux qe through the system is determined by 
the total precipitation inflow normalized by the dimensions 
of the flow layer 



Qprcc _ 



/ ^prcc \ 



(15) 



where Q pr0 c is the total precipitation inflow, g proc = 
Qprec/(L pTec Dy) is the specific precipitation flux incident 
through the horizontal plane, fi the relative precipitation 
fraction that goes into the downslope moisture flow, Z/ prC c 
the spatial extent of the high-precipitation source region in 
x above the flow layer, D y the spatial extent of the source 
region into the plane of the drawing Figure 3 in y, and d the 
thickness of the flow layer (see Section 5.1, Bear 1988). The 
results are independent of the scale D y , which divides out 
of the total precipitation inflow Q pr oc. 
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For purposes of estimating the downslope flux, we sup- 
pose that the total precipitation inflow Q prec can be divided 
systematically between components: ground moisture flow 
Q m , streamflow Q s f, evapotranspiration loss Q e t, and im- 
mediate run off Qro- For the dry 2000s 

Q P rec = Qm + Qsf + Qct + Qio, (16) 

and for the wet 1980s 

Qprcc = Qm + Qsf + Qct + Qro- (17) 

The Tularosa Creek streamflow overresponse to the precip- 
itation change suggests that there was less relative evapo- 
transpiration during the wet period. Total evaporation is 
determined by the exposed surface water, which should not 
change much with a small change in precipitation. Although 
the transpiration physics is complicated, natural mecha- 
nisms tend to compensate for decreased available moisture 
in the root zone so that plant transpiration may remain 
nearly constant with changing precipitation (Salvucci 2001; 
Teuling et al. 2006), thus we take Q et = Q c t- About 5% 
leaves as immediate flood runoff in 3-7 large storm events 
per year so Q ro = 0.05Q prC c and Q ro = 0.05Q prC c (estimated 
using New Mexico Water Resources Data, USGS yearly re- 
ports). The highest precipitation from 1987 was 33% higher 
than for the dry 2000s, Q prC c = 1.33Q pre c from the average 
line in Figure 2c, but the maximum relative Tularosa Creek 
streamflow was actually larger Qsf = 1.434Q s f from Figure 
2b. 

Though streamflow responds to the underground satura- 
tion in nonlinear ways as we describe in Section 8, suppos- 
ing for overall estimation purposes only that the streamflow 
stays approximately proportional to the downslope mois- 
ture flow for the relatively small changes as may be caused 
by precipitation variations in similar downhill saturation 
conditions, we have Q m + Q s f = 1.434(Q m + Q s f). Using 

Q ro = 0.05Q prcc = 0.05 • 1.33Q P rec = 0.0665Q prcc , we can 
eliminate all the quantities for the wetter period in (17), 
giving 

1.2635Qp re c = l-434(Q m + Q sf ) + Q ot . (18) 

Then substituting Q ro = 0.05Q prC c in (16) and eliminating 
Qm + Qsf with (18) gives the relative evapotranspiration 

Q ct = 0.228Q prcc , leaving Q m + Qsf = 0.722Q proc . (19) 

Under normal downhill saturated conditions in the area 
system, about 15% of the precipitation appears as stream- 
flow so Q s f = 0.15Q pre c, giving the ground moisture flow 
Q m = 0.572Q prec for the moisture fraction /x = Q m /Q P rcc = 
0.572. Under the reduced downhill saturation conditions 
in the La Luz / Fresnal Watershed after 2002, the relative 
streamflow was Q s f = 0.047Q prC c, which gives an increased 
relative moisture flow of Q m = 0.675Q prO c for fj, — 0.675. 

The precipitation influx for the La Luz / Fresnal Water- 
shed is q plC c = 25 in/yr = 1.74 mm/d asymptotically as a 
base level for 2006 from Figure 2c. We estimate the spatial 
extent for the high-precipitation source region as L proc ~ 
12 km as the flattening top of the watershed above 7000 ft 
(2134 m) elevation partially shown in Figure 1. A repre- 
sentative flow-layer thickness seems most uncertain, but we 
assume a typical spacing of the impervious limestone strata 
of d « 40 m characteristic of parts of the mountain system. 
With our guess and using the intermediate moisture fraction 
/i = 0.624, we obtain qe = 0.326 m/d from (15) reasonably 
consistent with the fluxes we obtained for the reference sat- 
uration degree S = 0.30. We must point out however that 
by adjusting the parameters within acceptable ranges, we 
obtain a range of fluxes and corresponding reference satura- 
tions. 



The downslope moisture flow speed is defined Ve = 
qe/(nS) from (12). Taking a porosity n = 0.17 for poorly 
sorted soils, with the flux qi = 0.326 m/d and intermediate 
saturation degree S = 0.6, we obtain the flow speed Ve = 3.2 
m/d. 

6. Physical Interpretation 

Steady equilibria like those derived in Section 2 repre- 
sent the long-term downslope profiles of saturation degree. 
Time series as represented by (1) and (3), or by the gen- 
eral Richards' Equation (4), or the ID equations (5) and (6) 
must tend to the nearest equilibrium with a characteristic 
timescale determined by the flow speeds for the problem. 

The ID steady lateral downslope equilibria are deter- 
mined by two free parameters: the flux qt and a boundary 
condition. Streamflows stay approximately uniform below 
the springs, but are absorbed into the ground within a few 
km beyond the foot of the mountains shown on the topo 
map in Figure 1, where the geology changes from sandy-clay 
limestone stratigraphy to alluvial fill. There the hydrology 
should change correspondingly from closed to open. In the 
basin, moisture pools with a variable deep water table and 
thus can exert no pressure back upslope. With an open 
boundary, no change in the downhill saturation is imposed 
and the downslope flow must follow the upper-hillslope form 
seen in all of the solution curves, a constant flow speed 
Vi — Ve and saturation degree S — S for the given flux 
qt from (12), like the solution reported in one study (Philip 
1991). With a change of precipitation inflow for a new flux 
qe, the flow speed Ve and saturation degree S must adjust to 
the new steady equilibrium defined in Section 4 with Figure 
7. 

At the onset of a new saturated zone, as may be pro- 
duced with a renewed sufficient continuous streamflow, the 
flux immediately upslope qe defined in (5) must be greatly 
decreased and may even become negative by the imposed 
dS/d£ positive singularity due to the large factor dip/dS near 
S = 1. The downslope flux is slowed or may even reverse its 
direction representing a net flux back upslope. The reduced 
downslope flux produces an increased upslope saturation ac- 
cording to the ID continuity condition (6). In effect the 
downslope moisture flow is blocked, and the blockage must 
propagate back upslope as moisture piles up behind it lead- 
ing to equilibrium solutions (9) with profiles like those shown 
in Figure 5, which represent a return to the equilibrium 
constant downslope-flux condition. We may say that mois- 
ture backs up behind the capillary head jump introduced in 
a zone of fixed saturation degree through a pooling height 
determined by the balance between the downslope gravity- 
driven flux out the bottom of the blockage and the influx 
into the top. If the saturated zone is not sufficiently main- 
tained, nonequilibrium or discontinuous behavior is possible 
and the flow may again return to the reference saturation 
below the saturated zone depending upon conditions further 
downslope. 

During the transition from normal streamflow to no 
streamflow or from a normally saturated downhill bound- 
ary condition to an open boundary condition, a temporary 
increase in the total moisture outflow into the basin below 
must occur. The normal flux out into the basin will return 
as the dryed system reaches a new equilibrium on the ad- 
justment timescale. The adjustment might be evidenced as 
a locally increased water table in the basin, except that the 
Tularosa Basin is so large that such changes may be diffi- 
cult to detect and probably could not be distinguished from 
other effects, like precipitation change or distributional inho- 
mogeneity. The Tularosa Basin surface area is about 13,500 
km compared to the La Luz / Fresnal Watershed, which is 
about 200 km 2 . Also it seems that adequate basin water- 
table measurements are not available. 
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7. Generalizations to 3D 

In real 3D conditions, moisture can travel out of the 
mountain source area along a number of independent paths, 
with different inflows, downslope variations in conditions, 
and differing downhill saturation zones on each path. The 
ID approximation should remain applicable in the idealized 
flow layer of Figure 3 even with changing material prop- 
erties downslope in i. As the equilibrium equations (5), 
(10), and (11), with qt constant in I are independent of the 
porosity n, a changing porosity downslope does not affect 
the downslope saturation solutions, but does affect the flow 
speed Vi(l) through (12). A changing hydraulic conductiv- 
ity or capillary pressure downslope in I does not affect the 
ID assumptions either since it only changes the material 
functions that enter (9), so solutions for S(i) should be like 
what is found in the uniform cases. 

A slowly changing downslope flux qt , as produced with an 
added precipitation influx, or evapotranspiration or flow- 
layer losses, in a thin flow layer of constant cross section 
must produce a slightly changed saturation degree downs- 
lope S(£) in the integrated (9), but should not alter the 
essential solution properties. Such variations should pro- 
duce systematic modifications to the pooling heights above 
saturated zones. 

With variations in conditions across the flow layer or a 
varying cross section or inclination downslope, non-lateral 
components to the flux vector may be introduced and the 
ID approximation break down, consistent with the belief 
that underground geological and topographical structure 
gives rise to deeper saturated structure like perched aquifers 
(Jackson 1992; Neuweiler and Cirpka 2005). In principle, 
the four equilibrium equations from (3) V • q = 0, and the 
three components of (1) define the three components of q 
and S throughout the 3D volume of the flow layer. The qt 
flux equation (5) still represents the downslope component 
in the general 3D conditions and must be preserved in in- 
tegrated form in the equilibrium solutions. Thus we argue 
that in 3D conditions where a nonnegligible downslope flux 
covering the cross section exists down the flow layer, the 
basic physical conservation principles from the ID solutions 
should still apply, and flow backup similarly above down- 
hill saturated zones that cover the cross section of a closed 
flow layer. Hurley and Pantelis [1985] define special forms 
of variations that still preserve the ID flow vector. 

The ID approximation breaks down too in the bound- 
ary zone if the saturating process by streamflow produces a 
non-uniform moisture distribution underground. Variations 
in the upslope pooling properties can occur if the under- 
ground is inhomogeneously partially saturated rather than 
being uniformly filled with a single saturation degree, as ide- 
alized in our ID solutions. While the two cases of a fully 
saturated boundary zone covering the cross section of the 
flow layer with maximum pooling above and an open bound- 
ary with no pooling are still correctly represented, partial 
coverage may change the properties of the downslope flow 
blockage and modify the pooling heights between the two 
extremes. 

8. Conclusions 

The hypothesis that water leaves the mountains mainly 
in a large-scale near-surface moisture flow is idealized by our 
ID Richards' Equation model for lateral downslope steady 
flow in a uniformly filled thin closed layer of constant cross 
section and inclination. The Richards' Equation represents 
general saturated/unsaturated underground moisture dy- 
namics. The ID steady flow equation is directly integrat- 
able and exhibits a constant downslope flow speed and mois- 
ture content determined by the constant downslope flux. At 
a boundary zone of increased saturation degree, the flow 



speed is slowed and moisture content increased back ups- 
lope through a pooling height, characteristized by the equi- 
librium condition that the flux, the product of the flow speed 
and moisture content, is constant. Surface streamflow must 
produce zones of increased saturation degree that affect the 
underground saturated boundary condition. Such a strong 
solution dependence on the form of the lower boundary con- 
dition suggests that mountain moisture content in hillslopes 
above deep alluvial basins may be strongly affected by lower 
streamflow. 

In our numerical integrations, the pooling height is more 
than 10 km vertically projected behind a fully saturated 
zone in all but the sandiest soil types, but decreases rapidly 
with a decreasing boundary saturation degree. The other 
cited studies of downslope moisture flow in a thin layer also 
show a constant flow speed and saturation degree or pool- 
ing at the bottom asymptotically in their time series, but 
our pooling heights appear to be much larger than what is 
found in the other studies. The large scale arises mainly be- 
cause material capillary head functions exhibit singular be- 
havior near full saturation, and all of the cited studies that 
show downstream pooling use different approximations for 
the material hydraulic conductivity or capillary-head func- 
tions to obtain a more tractable problem. 

Spring outflow must be a direct result of the coalescence 
of underground moisture in irregular topography near or up- 
slope from the spring (Freeze 1972; Fipps and Skaggs 1989). 
With a reduced downhill saturation, the saturation degree 
outside the streambed and back up to the spring may be 
greatly decreased from a near 100% pooling-determined sat- 
uration degree to a 30% reference saturation degree for rep- 
resentative conditions in the La Luz / Fresnal Watershed. 
Thus the drop of 69% in spring outflow, as the unaccounted 
loss seen at Alamogordo's upstream spring-box diversions, 
seems a reasonable possibility for the downhill drying, taking 
the spring outflow to be simply proportional to the satura- 
tion degree in the spring locale. Though individual spring 
response to changing conditions may be nonlinear, without 
knowing the specifics, the general suggestion from studies 
is that the central hillslope outflow increases with the lo- 
cal water-table height (cited references and Weyman 1973). 
The mountain saturation above may be significantly, but 
probably less compromised by downhill drying, as presum- 
ably not all of the moisture paths out of the mountain source 
area are affected. Also the adjustment time for the moun- 
tains above is longer and effects there should become evident 
only more slowly. 

The rate at which moisture moves downslope from a dis- 
appearing saturated zone or backs up behind a newly estab- 
lished one might be estimated supposing that an unmain- 
tained saturated zone will dissipate downslope at the av- 
erage flow speed, or flow backup behind a newly formed 
saturated zone at the same rate. The actual timescale for 
depletion of the La Luz / Fresnal Watershed seen in Figure 
2a and known from the history of the Alamogordo project 
was about 3 — 10 years, which for the flow speed Ve = 3.2 
m/d is the time required to travel 3.5 — 11.7 km, not out 
of range for the spring distances from the mountain-basin 
outlet shown in Figure 1. The estimate is consistent with 
the saturated-flow estimate described at the end of Section 
2. We note however the large uncertainty in the flow speed 
due especially to the uncertainties in the saturated hydraulic 
conductivity in real field conditions. 
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